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We propose a minimal nonlinear model of brittle crack propagation by considering 
only the motion of the crack-tip atom. The model captures many essential features of 
steady-state crack velocity and is in excellent quantitative agreement with many-body 
dynamical simulations. The model exhibits lattice-trapping. For loads just above this, the 
crack velocity rises sharply, reaching a limiting value well below that predicted by elastic 
continuum theory. We trace the origin of the low limiting velocity to the anharmonicity 
of the potential well experienced by the crack-tip atom. 
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Recent molecular-dynamics (MD) simulations of crack propagation [1][2], as well as 
experimental studies [3] [4], have reflected growing interest in the dynamical aspects of 
brittle fracture, including the approach to a steady (or quasi-steady) state, the buildup of 
coherent excitation near the crack tip [1][2], and the subsequent onset of instabilities [5] [6]. 
In all of these works, it is fair to say that a coherent, quantitative understanding of the 
limiting velocity dependence on the local field has not yet been advanced, though many 
good suggestions have been made [7] [2]. Here, we propose a minimal, one-atom, nonlinear 
model for describing brittle fracture, which we call the "Einstein Ice-Skater" (EIS) model. 

By closely observing movies of MD simulations of brittle crack propagation in a two- 
dimensional (2D) triangular lattice, under tensile (transverse, or Mode I) loading and 
at zero initial temperature, we noticed that cracks appear to advance as a sequence of 
essentially one-particle moves. Along the natural cleavage direction separating a pair of 
close-packed planes (lines in 2D), bond-breaking events are well separated in time [8], which 
can be characterized as a zigzag, ice-skating kind of motion between the two lines of atoms. 
When a bond breaks, the forward partner moves ahead, approximately along the former 
bond direction, while the rearward partner swings back to its final equilibrium position (see 
Figure 1). This led us to speculate that the steady state velocity of a brittle crack could 
be well approximated by a single-particle Einstein cell model, where the mobile crack-tip 
atom (the EIS in Figure 1) moves in a field of six immobile neighbors (the sixth, with whom 
the bond has just been broken, is assumed to be beyond the range of interaction). The 
bond-breaking event launches the EIS approximately along the bonding direction. This 
compressive, nonlinear event results in a shearing motion along the transverse pair of close- 
packed lines at ±60° to the propagation direction, and gives rise to the local vibrational 
excitations that build up around the crack tip and move coherently with it [1] [8] . 

For sufficiently large strains, the EIS reaches a point that stretches the next bond to 
breaking after a time tbreak since the last bond- breaking event. The pattern then repeats - 
to the other side of the ice-skating phase - and the crack has then advanced by one nearest 
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neighbor spacing r along the forward direction in the time 2tb rea k- The crack velocity is 
thus given by 



To find tbreak, we start from the configuration of the EIS and its five connected nearest 
neighbors and solve for the time dependence of the distance r i between the EIS and its 
neighbor #1; tbreak is the first time that roi reaches the breaking point r max . The equation 
of motion for the position r of the EIS (atomic mass m) is 



which can be solved given the pair potential 0(roi) = </>(l r o — r i|) an d the initial conditions. 
This equation is not trivial to solve, even for harmonic potentials, but can be solved 
numerically. We first assume the initial EIS coordinates x = a/2, y = and velocities 
x = y = (the initial velocities of steady-state crack-tip atoms in full MD simulations are 
observed to be indeed very small). With r a = 1, the six immobile neighbors are assumed 
to be located at (-a/2,1/2), (a/2,1), (3a/2,l/2), (3a/2, -1/2), (a/2 + a e, -1/2), and 
(—a/2 — aoe, —1/2), where ao = v^3/2, a = ao(l + e), and e is the uniaxial strain in the 
transverse direction to crack propagation. (See Figure 1.) 

We can obtain a crude estimate for tbreak by imagining that the EIS starts at the 
turning point of its motion in the final harmonic equilibrium well. The bulk Einstein model 
is characterized by a frequency of uje = a/3^o, where uo is the fundamental frequency given 
by tjilVq = 4>"{yq). Hence, if the time tbreak is one-half the period (from one turning point 
to the other at bond breaking), then w crac k = V3r uo/27i. Since the triangular-lattice 
shear-wave speed c s = ^/3/8r u)o (which is very close to the Rayleigh, or surface wave 
speed [9]), f C rack/c s ~ V^/tt = 0.45, independent of the anharmonicity of the potential. 
Since the effective frequency of a stretched anharmonic bond decreases (actually to zero 
at the inflection point), the crack velocity in the anharmonic case should be lower. 



"^crack = To/(2tbreak) • 



(1) 
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To go beyond this estimate, we investigated two kinds of attractive snapping-bond 
potentials: harmonic (HSB) and anharmonic (ASB), the latter based on the Morse poten- 
tial 

<f>(r) = [l-e- Q(r - 1) j 2 /(2a 2 ) . (3) 

Here we scale the distance by ro and the energy by mr qUq; a is the repulsive parameter (the 
familiar Lennard- Jones 6-12 potential is closely approximated by a = 6; most materials can 
be represented by 4 < a < 6). The ASB potential is obtained from Eq. (3) by flattening 
it out at r max = 1 — In (1 — y/x)/ a i* 1 the attractive region (beyond the minimum of the 
Morse potential); the cohesive energy is then x/(2ct 2 ), for X < 1- The ASB force jumps 
discontinuously at this point from a negative value to zero - hence the term "snapping- 
bond." For small displacements about r = 1, Eq. (3) is approximately harmonic, (r— 1) 2 /2. 
The HSB potential cuts off at the same energy as the ASB, but at r^ ax = 1 + < r max . 

We find that the range and maximum attractive force of the potential are the essential 
parameters that govern the crack velocity. 

Our choice of snapping-bond potentials makes precise the definition of the distance 
beyond which a bond is considered "broken," an ambiguous concept for completely con- 
tinuous potentials. Since our goal is to compare this EIS model with a fully dynamical 
system, a well-defined breaking point for both is a distinct advantage. The fully dynam- 
ical systems we compare with are rather restrictive, namely, close-packed lines of atoms 
of width w = 4,8,16, and 64, with the outer two clamped, and the inner free to move; 
moreover, only nearest-neighbor interactions are considered. (Strips were typically 200r o 
in length; steady-state propagation is attained well within 10% of that length.) 

For this thin-strip, fixed-grip geometry, the critical Griffith strain €q for initiating 
forward crack motion can be computed by equating the potential energy in two transverse 
sections of the strip of height ro/2: one far behind the crack with all bonds in equilibrium, 
except for the one broken bond, and the other far in front, with all bonds equally stretched. 
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The Griffith criterion is obtained from 



(w - l)0(n) = 0(r 2 ) = X /(2a 2 ) , (4) 

where r\ is the elastically stretched bond [r\ = a 2 + 1/4) and is the broken bond across 
the gap of the relaxed crack. The Griffith criterion ec is thus 




An intriguing aspect of the EIS model is the straightforward emergence of the lattice 
trapping phenomenon [10]: unless the strain exceeds a value well above ea, the distance 
between the EIS and its neighbor #1 will not reach r max . The strain must therefore exceed 
€q by a barrier amount of overstrain that is a characteristic of the atomistic nature of the 
crack tip, and which can only be evaluated atomistically. In Figure 2, we show our results 
for the crack-tip velocity (in units of c s ), as a function of the strain, for the EIS model and 
for the fully dynamical w = 4 strip (a = 6, x = 1/2)- The EIS model agrees to within 
10% of the velocity with the MD results - remarkable for such a simple model. 

However, the lattice-trapping strain is underestimated by 13% for the anharmonic 
and 12% for the harmonic system, which is most clearly due to neglected correlations 
with farther neighbors in the EIS model. For the anharmonic system, the onset of crack 
motion for the fully dynamical w = 4 strip occurs at a crack velocity of about 30% of the 
shear-wave speed, while for the harmonic system, the crack starts at about 50% of the 
shear-wave speed. Under further loading, the crack-tip velocity increases roughly linearly 
with strain but with a higher slope for the harmonic than for the anharmonic system. 

To compare our EIS results to MD simulations and experiments we rescaled the wider 
system strains by the Griffith strain (e^ ~ and found good agreement, except for 

slight, but systematic increases in the lattice trapping strain with size for harmonic poten- 
tials. We can understand this by noting that wide anharmonic systems, where stretched 
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bonds weaken, are more compliant and tend to have local strains near the crack tip that 
are closer to those in the narrow strips. On the other hand, harmonic bonds do not weaken 
with stretching, so that the global strain is spread more democratically across the system. 
We emphasize that, even in wide systems where the global strain can be arbitrarily small, 
the fact that local strains near the crack tip are large (of order 10%, as in the narrow-strip 
case) is a significant reason for the success of the EIS model. 

We find that crack velocities in anharmonic systems are essentially independent of 
the anharmonicity parameter, at least over the range 4 < a < 6; in fact, the curves for 
a = 4 and 5 practically overlap. As the cohesive strength \ decreases from 1/2 down 
to 1/8 (along with the range of the potential), crack velocities in anharmonic systems 
show a slight increase (~ 10%) in ultimate slope and greater variability in the jump-off 
lattice-trapping strain. (In the limit x ~~ 0, of course, the harmonic limit is approached 
[7].) In general, velocities in anharmonic systems are lower than in harmonic ones, show 
less variation with strain, and exhibit relatively lower lattice trapping (when the strain is 
scaled by e^). Similar trends are exhibited in the full MD simulations, including those 
using full, continuous (rather than discontinuous snapping-bond) potentials [1] [2] [7] [8] , 
and those for systems much wider than w = 4. Again, the principal differences are in the 
lattice-trapping strains. We can therefore conclude that the EIS approximates very well 
the crack-tip atomic motion, just as our intuition from larger-scale MD simulations had 
suggested. 

Our minimal EIS model indeed confirms speculations about the correlation of the 
limiting steady-state crack-tip velocity and anharmonicity [2] [7] [8]. The more "realistic" 
anharmonic interactions give steady-state crack-tip velocities that never exceed 0.4 of the 
Rayleigh speed, in excellent agreement with experimental observations [3] [4]. With the 
EIS model, the origin of this low speed can clearly be attributed to the smaller attractive 
force on the crack-tip atom at the point of bond breaking, as compared to the harmonic 
(or linear elastic) analysis. 
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Under loading, the thin-strip MD crack-tip velocity in Figure 2 jumps sharply at the 
lattice-trapping strain to a slowly rising linear regime, and then once again rises sharply 
at a strain of 0.15 = 1.3ec. Close examination of atomic configurations revealed that this 
second rise is associated with two instabilities: the first is a wake of large-amplitude surface 
(Rayleigh) waves behind the crack tip; at somewhat higher strains, the crack begins to 
jump from the central channel to one of the side channels next to the fixed-grip atoms 
(see Figure 1). We have observed this zigzag propagation by plus or minus one channel 
in much wider systems, where, at even higher strains, dislocations are emitted, followed 
immediately by branching. Dynamical instabilities such as these divert energy from brittle 
bond breaking, causing the crack-tip velocity to drop rather than rise. Dislocation emission 
and real crack branching are, of course, forbidden processes in the artificially narrow 4- wide 
strip, and are completely absent in the one-particle EIS model. 

Finally, the hysteresis under unloading and healing up of the crack can be obtained 
from the EIS. To do this, we simply detect when the 6-neighbor model reconnects the bond 
between the EIS atom and neighbor #6, rather than opening up the crack in the forward 
direction. This occurs soon below eo for the anharmonic potential (namely, 0.98cg), but 
substantially lower for the harmonic potential (0.85ec). Crack propagation and crack 
healing are thus quite asymmetric processes. 

In conclusion, the Einstein Ice-Skater model of brittle crack propagation is able to pre- 
dict quantitatively the steady-state crack velocity under loading, including lattice trapping, 
as well as hysteresis upon unloading and crack healing. The maximum velocity achieved 
in full MD simulations as a function of strain is principally limited by the anharmonicity 
in the attractive region of the pair potential, which is captured by the EIS; however, it 
is also affected by instabilities that involve collective motion (energy buildup, dislocation 
emission, and branching), which is inaccessible to the one-particle EIS model. Neverthe- 
less, this simple EIS model allows us to explain, in quite satisfactory quantitative fashion, 
the effect of nonlinear motion of the crack-tip atom on the limiting crack velocity. 
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FIGURE CAPTIONS 

1. Initial atomic coordinates for crack propagation in a triangular-lattice strip, four close- 
packed lines wide; the outer two lines of atoms are fixed, while the inner two are 
mobile. Heavy lines indicate equilibrium (nearest-neighbor) bonds of length r$ = 1; 
heavy dashed lines are slightly stretched, nearly vertical bonds; light lines are bonds 
elastically stretched to length r\ ~ 1 + 3e/4 by the uniaxial strain e in the x-direction; 
the light dashed line is a just-broken bond with neighbor #6. The EIS atom is 
indicated by the large open circle: it moves initially approximately in the direction of 
the arrow, stretching the bond with neighbor #1 until breakage, then heads toward 
its final equilibrium position (small circle). 

2. Crack velocity (in units of shear- wave speed c s ) as a function of strain for the anhar- 
monic snapping-bond (ASB) and harmonic snapping-bond (HSB) potentials. Results 
for the EIS model are shown for Morse parameter a = 6 and cohesive bond-strength 
X = 1/2, along with w = 4 strip MD simulations (closed circles for ASB and open for 
HSB). 
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